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Abstract 

Mathematical  modeling  plays  an  important  role  in  fuel  cell  design.  A  comprehensive  review  of  the  mathematical  modeling  of  proton  exchange 
membrane  fuel  cells  is  first  conducted.  It  is  found  that  the  results  computed  by  different  models  in  the  literature  often  agree  well  with  the  experimental 
data.  This  stimulates  the  present  authors  to  carry  out  a  comprehensive  parameter  sensitivity  examination.  In  this  first  paper  a  three-dimensional, 
two-phase  and  non-isothermal  model  is  developed,  and  numerical  simulations  for  a  basic  case  is  performed,  the  results  of  which  are  regarded  as 
the  reference  for  further  sensitivity  examination.  All  the  parameters  needed  for  the  simulation  are  provided  in  detail.  In  the  companion  paper  (Part 
II),  the  results  of  the  parameter  sensitivity  analyses  and  discussion  of  model  validation  are  provided  in  detail. 

©  2006  Elsevier  B.Y.  All  rights  reserved. 
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1.  Introduction 

The  proton  exchange  membrane  fuel  cell  (PEMFC)  is  consid¬ 
ered  to  be  a  promising  power  source,  especially  for  transporta¬ 
tion  and  stationary  cogeneration  applications  due  to  its  high  effi¬ 
ciency,  low-temperature  operation,  high  power  density,  fast  start¬ 
up,  and  system  robustness.  In  the  last  decade  a  great  number  of 
researches  have  been  conducted  to  improve  the  performance  of 
the  PEMFC,  so  that  it  can  reach  a  significant  market  penetration. 
Ref.  [1-5]  are  the  examples  of  very  recent  publications.  In  this 
regard,  an  optimization  study  of  the  PEMFC  plays  an  important 
role.  One  of  the  important  tools  in  the  optimization  study  of  fuel 
cell  performance  is  computational  modeling,  which  can  be  used 
to  reveal  the  fundamental  phenomena  taking  place  in  the  fuel  cell 
system,  predict  fuel  cell  performance  under  different  operating 
conditions,  reveal  the  distribution  details  of  various  dependent 
variables  and  optimize  the  design  of  a  fuel  cell  system  [6] . 
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The  processes  in  the  fuel  cells  are  very  complicated  because 
of  the  very  tight  coupling  between  electrochemical  and  transport 
processes.  For  modeling  of  a  single  fuel  cell,  the  parameters  of 
electrochemical  kinetics,  fluid  flow,  mass  transfer,  heat  transfer 
and  species  transfer  are  included.  Because  the  number  of  the 
flow  channels  in  the  bipolar  plate  is  quite  large,  and  there  are 
seven  functional  regions  across  the  fuel  cell,  with  the  present- 
day  models  available  to  most  researchers  are  often  impossible 
to  use  to  numerically  simulate  the  whole  fuel  cell.  Therefore,  a 
typical  element  is  usually  separated  from  the  whole  fuel  cell  in 
the  computational  domain.  Although  this  element  only  covers 
part  of  the  fuel  cell,  it  includes  all  functional  parts  of  the  fuel 
cell:  from  the  anode  channel  to  the  cathode  channel.  Thus,  based 
on  the  assumption  that  the  process  is  periodic  from  channel  to 
channel,  such  an  element  may  be  regarded  as  the  representation 
of  the  entire  fuel  cell.  This  kind  of  numerical  simulation  may  be 
called  a  typical  unit  simulation  of  a  single  fuel  cell.  The  work 
reviewed  and  presented  in  this  paper  belongs  to  this  category  of 
simulation. 

Because  of  the  complexity  of  the  process,  there  are  more 
than  ten  empirical  or  experimental  parameters  involved  in  the 
physical  modeling  of  the  PEMFC.  To  validate  the  physical  and 
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Nomenclature 

a  water  activity 

A  area  (m2) 

As  specific  area  of  catalyst  layer  (m-1) 

c  molar  concentration  (mol  m-3) 

D  diffusion  coefficient  (m2  s-1) 

F  Faraday’s  constant  (C  mol- 1 ) 

hm  evaporation  and  condensation  rate 

H  height  (m) 

i  reaction  rate  (A  m-  3 ) 

/re f  reference  exchange  current  density  (A  m-2) 

I  current  density  (A  m-2) 

k  thermal  conductivity  (W  m-1  K-1) 

K  electrode  absolute  permeability  (m2) 

L  length  (m) 

M  molar  mass  (kg  mol  - 1 ) 

n  electron  number  of  electrochemical  reaction 

p  pressure  (Pa) 

R  universal  gas  constant  (8.314  J  mol-1  K-1) 

.v  liquid  water  saturation 

S  source  term  of  governing  equations 

T  temperature  (K) 

u  velocity  vector  (ms-1) 

V  potential  (V) 

w  velocity  at  z-direction  (ms-1) 

W  width  (m) 

v,  y,  z  coordinate  (m) 

X  species  mass  fraction  (dimensionless) 


Greek  symbols 
a  transfer  coefficient 

s  porosity 

f  stoichiometric  flow  ratio 

v)  overpotential  (V) 

k  electrical  conductivity  (S  m- 1 ) 

X  membrane  water  content 

/j  viscosity  (kg  m- 1  s- 1 ) 

v  kinetic  viscosity  (m2  s-1) 

p  density  (kg  m- 3 ) 

o  surface  tension  (N  m- 1 ) 

0  potential  (V) 

co  species  molar  fraction 


Subscripts  and  superscripts 
a  anode 

av  average  value 

c  cathode 

cc  land  area 

eh  channel 

ct  catalyst  layer 

d  diffusion  layer 

eff  effective 

g  gas 

h  hydrogen 

in  inlet 


k  species 

1  liquid 

m  membrane 

o  oxygen 

oc  open  circuit 

r  relative  values 

ref  reference  values 

s  solid;  specific 

sat  saturation 

tot  total 

w  water 


numerical  models,  comparison  with  some  experimental  data  is 
highly  desirable.  For  the  fuel  cell  performance  description,  the 
polarization  curve,  or  voltage-current  curve,  is  one  of  the  most 
important  final  outcomes  of  numerical  simulation.  For  a  com¬ 
parison  with  experimental  results,  the  most  widely  cited  test 
data  of  the  polarization  curve  in  previous  literature  are  the  ones 
presented  by  Ticianelli  et  al.  in  [7,8].  By  careful  review  of  the 
existing  literature  available  to  the  present  authors,  it  was  found 
that  different  models  with  different  parameters  were  adopted  in 
simulations,  while  the  final  outcome,  i.e.,  the  curves  for  the  fuel 
cell  voltage  versus  current  were  often  almost  the  same.  What 
was  more  interesting  was  the  fact  that  the  comparison  of  dif¬ 
ferent  simulation  results  with  the  same  test  data  [7,8]  often 
showed  good  agreement.  This  situation  stimulates  the  present 
authors  to  conduct  a  parameter  sensitivity  examination  for  two 
purposes:  Firstly,  to  reveal  what  parameters  have  the  most  signif¬ 
icant  effects  on  the  V-I  curve;  secondly,  to  examine  whether  the 
comparison  with  test  data  of  the  V-I  curve  is  enough  to  validate 
a  model. 

This  paper  is  the  first  part  of  a  two-part  study  on  the  param¬ 
eter  sensitivity  study  and  a  discussion  of  the  validation  model. 
The  rest  of  the  paper  is  organized  as  follows.  In  order  to  have 
a  clear  understanding  of  the  above-mentioned  situation,  a  com¬ 
prehensive  review  of  numerical  models  will  first  be  conducted. 
Then,  a  three-dimensional  two-phase  and  non-isothermal  model 
will  be  presented  for  the  parameter  sensitivity  study.  Numerical 
methods,  a  solution  flow-chart  and  grid-independence  examina¬ 
tion  for  the  proposed  model  will  be  presented.  The  numerical 
simulation  results  for  a  basic  case  will  be  shown  which  as  the  ref¬ 
erence  for  a  further  parameter  sensitivity  examination.  Finally 
some  conclusions  will  be  made.  The  results  of  the  parameter  sen¬ 
sitivity  study  and  a  detailed  discussion  of  the  model  validation 
will  be  presented  in  the  companion  paper. 

2.  Review  of  existing  PEM  fuel  cell  models 

The  schematic  of  a  PEMFC  is  shown  in  Fig.  1.  The  anode 
gas  channel  (in  a  bipolar  plate),  the  anode  diffusion  layer,  the 
anode  catalyst  layer,  the  ion-conducting  membrane,  the  cath¬ 
ode  catalyst  layer,  the  cathode  diffusion  layer,  and  the  cathode 
gas  channel  are  components  of  the  PEMFC.  Bipolar  plates  act 
as  electron  collectors.  Anode  gas  channels  supply  the  fuel  cell 
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Fig.  1.  Schematic  view  of  a  PEMFC. 

with  reactants,  which  are  transported  by  convection  and  dif¬ 
fusion  throughout  the  anode  gas  channel.  The  electronically 
conducting  porous  diffusion  layers  allow  for  more  or  less  even 
distribution  of  the  reactants  over  the  anode  and  cathode.  Hydro¬ 
gen  oxidation  and  oxygen  reduction  reactions  are  considered  to 
occur  only  within  the  anode  catalyst  layers  and  cathode  catalyst 
layers,  respectively.  The  catalyst  layers  also  provide  channels 
for  transfer  of  electrons  and  protons.  The  ion-conducting  mem¬ 
brane  spatially  separates  the  fuel  cell  into  two  parts  and  only 
allows  the  transport  of  water  and  protons. 

2.7.  Classification  of  the  PEMFC  simulation  models 

A  number  of  simulation  models  have  been  developed.  These 
models  can  generally  be  characterized  by  the  computational 
scope  of  the  model.  One  kind  of  model  focuses  on  a  specific 
part  or  parts  of  the  fuel  cells  [9],  such  as  the  cathode  catalyst 
layer  [10-14],  the  electrode  [15-18],  the  gas  diffusion  layer 
[19],  the  membrane  electrode  assembly  [20,21],  or  the  ion¬ 
conducting  membrane  [22-24].  Ref.  [25]  reviewed  the  Nation 
membrane  structure  and  properties.  These  models  are  useful. 
However,  they  cannot  provide  a  complete  picture  of  the  fuel 
cell.  The  other  kinds  of  models  include  all  parts  of  a  fuel  cell, 
from  one-dimensional  and  single-phase  to  three-dimensional 
and  two-phase. 

The  most  prominent  earlier  work  was  from  Bernardi  and  Ver- 
brugge  [21,26,27]  and  Springer  et  al.  [20,28],  who  developed 
one-dimensional  models.  Later,  Eikerling  et  al.  [29],  Baschuk 
and  Li  [30],  Rowe  and  Li  [31]  and  Maggio  et  al.  [32]  devel¬ 
oped  other  one-dimensional  models.  These  models  can  predict 
the  cell  V-I  performance  in  the  low  and  intermediate  current 
density  ranges  with  reasonably  good  agreement  of  the  test  data, 
but  fail  to  reproduce  the  concentration  polarization  in  the  polar¬ 
ization  region  [33].  Therefore,  two-dimensional  models  were 
presented  by:  Kulikovsky  et  al.  [1 1],  Nguyen  and  White  [34],  Yi 
and  Nguyen  [35],  Gurau  et  al.  [36],  Kazim  et  al.  [37],  Singh  et 
al.  [38],  He  et  al.  [39],  Wang  et  al.  [40],  Natarajan  and  Nguyen 
[41]  and  Hsing  and  Futerko  [42].  A  relatively  simpler  approach 
to  model  the  plane  formed  by  the  direction  across  the  fuel  cell 
and  the  direction  along  the  flow  channel  is  called  a  quasi-two- 


dimensional  model  [43,44].  In  these  models,  a  one-dimensional 
model  for  one  direction  is  coupled  with  another  one-dimensional 
model  for  another  direction  normal  to  the  previous  one  to  sim¬ 
ulate  the  profiles  in  the  plane.  The  two-dimensional  models  can 
only  simulate  the  plane  perpendicular  to  the  flow  channels  (in 
z-y  plane  of  Fig.  1)  or  the  plane  formed  by  the  direction  across 
the  fuel  cell  and  the  direction  along  the  flow  channel  (z-x  plane 
of  Fig.  1).  Therefore,  these  models  cannot  give  a  full  picture  of 
the  variations  in  the  temperature  and  reactants  in  a  typical  three- 
dimensional  element.  In  order  to  have  a  better  understanding 
of  how  the  actual  fuel  cell  performs,  it  is  necessary  to  have  a 
three-dimensional  model.  Three-dimensional  models  are  being 
developed.  Dutta  et  al.  [45,46],  Shimpalee  and  Dutta  [47],  Bem- 
ing  et  al.  [48],  Berning  and  Djilali  [49],  Jen  et  al.  [50],  Fi  et  al. 
[51],  Hu  et  al.  [52],  Um  and  Wang  [53],  Nguyen  et  al.  [54], 
Zhou  and  Fiu  [55]  and  Ju  et  al.  [56]  have  all  presented  three- 
dimensional  models.  Fike  the  quasi-two-dimensional  model, 
some  quasi-three-dimensional  models  are  developed  in  [57,58]. 

2.2.  Validation  issues  for  P  EM  fuel  cell  models 

As  indicated  above,  modeling  research  on  the  PEM  fuel  cell  is 
being  developed  quickly.  As  indicated  in  [59]  all  models  should 
be  validated  by  experimental  data  such  as  presented  in  [7,8] 
or  by  other  successful  models  [60].  From  the  literature,  most 
of  the  models  are  validated  by  experimental  data.  The  general 
method  used  is  the  polarization  curve,  i.e.,  V-I  curve  computed 
by  the  model  is  compared  with  the  polarization  curve  measured 
by  experiment.  If  two  V-I  curves  are  in  good  agreement,  the 
model  is  usually  considered  reliable.  There  are  two  methods 
of  comparison  with  test  data.  In  one,  the  authors  who  devel¬ 
oped  the  model  measured  the  experimental  data  themselves 
[8,31,56,61-66]  or  together  with  other  research  groups  [67]. 
However,  most  authors  take  the  other  route:  validating  their  mod¬ 
els  by  comparing  their  numerical  results  with  the  experimental 
data  published  by  other  research  groups  [10,30,52,60,67-79]. 
Our  discussion  is  focused  on  the  second  method.  We  have  found 
that  often  the  results  computed  by  different  models  agree  surpris¬ 
ingly  well  with  the  same  experimental  data.  It  should  be  noted 
that  many  empirical  parameters  involved  in  different  models  are 
often  very  different  or  even  rather  different  from  the  experimen¬ 
tal  data  used  for  validation  of  the  model.  For  example,  we  find 
that  there  are  at  least  fourteen  papers  which  used  the  experi¬ 
mental  data  published  by  Ticianelli  et  al.  [7,8]  to  validate  their 
models.  These  include  following  papers:  the  one-dimensional 
single-phase  model  of  [26,60,80],  two-dimensional  single-phase 
models  of  [35,37,81-83],  three-dimensional  single-phase  mod¬ 
els  of  [47,49,8 1 ,82,84]  and  three-dimensional  two-phase  models 
of  [85-87].  We  find  that  [7,8]  did  not  supply  all  of  the  parameters 
the  models  need.  These  two  papers  supplied  many  polarization 
curves  under  different  conditions.  Many  papers  used  the  same 
curve  under  conditions  corresponding  to  the  following  case: 
polymer-electrolyte,  20  wt.%  Pt,  50-nm  sputtered  Pt,  5  atm  cath¬ 
ode  pressure,  3  atm  anode  pressure  and  a  353  K  cell  temperature. 
Some  representative  parameters  used  in  ten  papers  are  listed  in 
Table  1,  which  shows  that  most  of  the  parameters  used  in  the  ten 
papers  are  different  from  each  other,  but  the  V-I  curves  of  the 
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Table  1 


Parameters  used  in  10  papers 


Unit 

[26] 

[35] 

[81,82] 

[83] 

[48] 

[84] 

[85,86] 

[87] 

L 

m 

0.0762 

0.07112 

0.07112 

0.0767 

0.01067 

0.0711 

0.05 

W 

m 

7.62  x  10“4 

1.0  X  10“3 

1.0  X  10“3 

Wcc 

m 

1.0  x  10“3 

1.0  x  10“3 

#ch 

m 

7.62  x  10“4 

7.62  x  10“4 

7.62  x  10“4 

5.0  x  10“3 

7.62  x  10“4 

7.62  x  10“4 

1.0  x  10“3 

Hd 

m 

2.6  x  1CT4 

2.54  x  10“4 

2.54  x  10“4 

2.54  x  10“4 

2.54  x  10“4 

2.54  x  10“4 

2.54  x  10“4 

4.16  x  10“4 

Ha 

m 

1.0  x  10“6 

2.87  x  10“5 

2.87  x  10“5 

2.87  x  10“5 

2.87  x  10“5 

2.87  x  10“5 

2.87  x  10“5 

2.01  x  10“5 

Hm 

m 

2.3  x  10“4 

2.3  x  10“4 

2.3  x  10“4 

2.3  x  10“4 

2.3  x  10“4 

2.3  x  10“4 

1.75  x  10“4 

1.25  x  10“4 

Tfi.rcf 

2  -1 
m  s 

6.5  x  10“5 

5.22  x  10“6 

5.22  x  10“6 

5.22  x  10“6 

^h,rcf 

2  -1 
m  s 

3.7  x  10“5 

2.63  x  10“6 

3.76  x  10“6 

3.76  x  10“6 

K 

2 

m 

1.76  x  10“H 

1.76  x  10”11 

1.76  x  10“n 

1.76  x  10“u 

1.76  x  10“u 

1.76  x  10”11 

1.76  x  10“n 

Ks 

S  m~ 1 

120 

53 

53 

Km 

Sm"1 

17 

[20]* 

[20]* 

[20]* 

[20]* 

17 

[20]* 

4  s  Crcf 

Am-3 

9.23  x  108 

5.0  x  10s 

5.0  x  10s 

5.2  x  108 

8.5  x  10s 

YCrcf 

Am~3 

5.0  x  102 

1.05  x  106 

1.0  x  102 

1.0  x  102 

1.1  x  102 

45 

£d 

0.4 

0.4 

0.4 

0.4 

0.4 

0.4 

0.3 

ota 

1.0 

0.25 

1.0 

1.0 

0.25 

1.0 

1.0 

ac 

0.5 

0.5 

0.625 

0.5 

0.5 

0.625 

0.5 

0.5 

Ch,ref 

mol  m-3 

56.4 

co,ref 

mol  m-3 

4.62 

3.39 

The  value  of  /cm  in  the  model  is  computed  by  the  method  in  [20]. 


ten  papers  all  agree  well  with  the  same  experimental  data.  This 
very  interesting  situation  stimulates  the  authors  to  conduct  the 
present  study.  In  the  following,  a  review  will  first  be  given  on 
the  previous  parameter  effects  study  and  the  model  validation 
issue  for  the  PEMFC. 

2.3.  Previous  parameter  effects  studies  and  model 
validation  discussion 

Many  researchers  have  made  an  investigation  of  the  effects 
of  parameters  on  the  performance  of  fuel  cells.  The  parameters 
which  have  been  investigated  include:  operating  parameters  (i.e. 
temperature,  humidity,  pressure,  flow  rate,  etc.),  design  param¬ 
eters  (i.e.  geometrical  parameters,  etc.),  physical  parameters 
(i.e.  porous  media  porosity,  permeability,  etc.),  electrochemical 
parameters  (i.e.  specific  area  of  catalyst  layer  timed  by  reference 
electrical  density)  and  other  parameters  (i.e.  CO  poisoning,  etc.). 

Berning  et  al.  [48]  investigated  the  effects  of  various  operat¬ 
ing  parameters  such  as  temperature,  pressure  and  stoichiometric 
flow  ratio  on  the  fuel  cell  performance.  In  addition,  geometri¬ 
cal  and  material  parameters  such  as  the  gas  diffusion  electrode 
thickness  and  porosity  as  well  as  the  ratio  between  the  channel 
width  and  the  surface  area  were  investigated.  Yi  and  Nguyen 
[88]  investigated  the  effects  of  the  gas  hydrodynamics  on  the 
performance  of  the  air  cathode  of  a  PEM  fuel  cell  in  contact 
with  an  inter- digitated  gas  distributor.  In  addition  the  effects  of 
pressure  drop  between  the  inlet  and  outlet  channels  of  an  inter- 
digitated  gas  distributor,  the  electrode  height,  and  shoulder  width 
on  the  average  current  density  were  investigated.  Pasaogullari 
and  Wang  [89]  developed  a  model  to  explore  the  two-phase 
flow  physics  in  the  cathode  gas  diffusion  layer.  The  simulations 
revealed  that  flooding  of  the  porous  cathode  reduced  the  rate 
of  oxygen  transport  to  the  cathode  catalyst  layer.  Furthermore, 
they  indicated  that  the  humidification  level  and  the  flow  rate  of 


reactant  streams  are  key  parameters  controlling  PEMFC  perfor¬ 
mance  and  two-phase  flow  and  transport  characteristics.  Jaouen 
et  al.  [16]  developed  a  one-dimensional,  steady-state  agglom¬ 
erate  model  to  study  the  nature  of  mass  transport  limitations 
in  the  PEMFC  cathode.  The  effect  of  the  active  layer  thick¬ 
ness,  oxygen  concentration  and  relative  humidity  of  the  oxygen 
stream  were  investigated.  Chen  et  al.  [90]  presented  a  two- 
dimensional,  along-the-channel  model  to  design  fuel  channels 
for  proton  exchange  membrane  (PEM)  fuel  cells.  The  analysis 
was  made  of  the  effects  of  some  operation  and  design  parame¬ 
ters,  such  as  inlet  velocity,  inlet  pressure,  catalyst  activity,  height 
of  channel,  and  porosity  of  gas-diffusion  layer.  Chu  et  al.  [91] 
and  Jeng  et  al.  [92]  made  an  investigation  of  the  effects  of  the 
change  of  the  gas  diffuser  layer  porosity  on  the  performance  of 
a  proton  exchange  membrane  fuel  cell.  Lee  et  al.  [93]  presented 
a  simulation  of  the  fluid  in  the  gas  channel  and  the  diffusion 
layer  for  the  effects  on  the  electrode  variables:  gas  diffusion 
layer  thickness,  porosity,  and  distribution  of  pore  size.  Lum  and 
McGurik  [94]  developed  a  model  of  the  cathode  of  a  PEMFC 
with  an  inter-digitated  gas  distributor  with  the  intention  of  study¬ 
ing  the  effects  of  various  operating  parameters  such  as  electrode 
permeability,  thickness  and  shoulder  width.  Kazim  et  al.  [95] 
developed  a  two-dimensional  mathematical  model  to  investigate 
the  effects  of  parameters  such  as  cathode  porosity,  inlet  oxygen 
mole  fraction,  operating  temperature  and  pressure.  Hwang  et  al. 
[96]  developed  a  three-dimensional  numerical  model  to  simulate 
the  transport  phenomena  on  the  cathodic  side  of  a  PEMFC  and 
compared  the  polarization  curves  of  the  inter-digitated  flow  field 
and  parallel  flow  field.  Recently,  Li  and  Sabir  [97]  presented  a 
review  of  the  state-of-the-art  for  different  bipolar  plates  in  PEM 
fuel  cells.  Meng  and  Wang  [98]  investigated  effects  of  electron 
transport  through  the  gas  diffusion  layer  in  detail.  Du  et  al.  [99] 
investigated  the  effective  proton  and  electronic  conductivity  of 
the  catalyst  layers.  Chan  and  Tun  [100]  investigated  the  effects 
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of  the  cathode  reference  exchange  current  density  multiplied  by 
the  area,  reference  oxygen  concentration  and  oxygen  diffusivity 
on  the  performance.  Sun  et  al.  [101]  applied  a  two-dimensional 
cross-the-channel  model  to  investigate  the  influence  of  the  gas 
diffusion  layer  property  and  flow-field  geometry  such  as  dif¬ 
fusion  layer  diffusivity,  diffusion  layer  conductivity,  channel 
width-to-area  ratio  and  diffusion  layer  thickness  on  the  local 
reaction  rate  in  the  PEMFC  cathode  catalyst  layer.  They  found 
that  when  the  PEMFC  uses  reformate  hydrogen,  the  PEMFC 
performance  drops  dramatically  due  to  CO  poisoning  as  the 
anode  flow  rate  increases.  Zhou  and  Liu  [55],  Chan  et  al.  [76], 
Baschuk  et  al.  [78,102],  Camara  et  al.  [103],  Zhang  et  al.  [104] 
and  Wagner  and  Giilzow  [105]  all  investigated  the  effects  of  CO 
poisoning  on  PEMFC  performance.  More  recently,  researchers 
started  to  focus  on  the  investigation  of  the  air-breathing  PEMFC 
models  [106-109]. 

In  the  above  references  the  major  focus  mainly  was  con¬ 
centrated  on  the  effect  of  some  individual  parameters,  the 
comparison  of  effects  from  different  parameters  and  the  sen¬ 
sitivity  of  the  parameter  variation  on  the  final  outcome  of  V-I 
were  not  their  main  purpose.  Here  the  parametric  sensitivity 
refers  to  what  degree  a  parameter  affects  the  PEMFC  perfor¬ 
mance,  especially  the  V-I  curve.  Recently,  some  researchers 
have  been  aware  of  the  importance  of  the  sensitivity  issue  and 
several  papers  were  published.  Stockie  et  al.  [110]  performed 
a  sensitivity  study.  It  was  revealed  that  some  geometrical  and 
operational  parameters  are  critical  to  fuel  cell  performance. 
Grujicic  et  al.  [Ill]  performed  a  sensitivity  analysis  to  deter¬ 
mine  the  effects  of  six  parameters,  including  thickness  of  the 
active  layer,  molar  diffusion  volume  of  oxygen,  molar  diffusion 
volume  of  water,  molar  diffusion  volume  of  nitrogen,  and  poros¬ 
ity  of  cathode,  on  the  performance  of  PEM  fuel  cell  based  in 
a  steady-state  single-phase  three-dimensional  electro-chemical 
model.  The  results  showed  that  the  performance  of  a  common 
PEMFC  is  strongly  affected  by  these  parameters.  But  the  per¬ 
formance  of  the  fuel  cell  specially  designed  by  their  optimizing 
method  is  essentially  unaffected  by  these  parameters.  Correa 
et  al.  [112]  presented  a  sensitivity  analysis  on  the  PEMFC 
stack.  They  classified  the  parameters  according  to  their  influ¬ 
ence  in  the  fuel  cell  stack  as:  insensitive,  sensitive,  and  highly 
sensitive. 

From  Table  1  and  the  above-referenced  papers,  it  can  be  found 
that  the  transport  and  electrochemical  parameters  involved  in  the 
numerical  modeling  have  a  wide  range  of  variation.  The  present 
authors  then  examined  how  the  values  of  parameters  are  obtained 
in  the  numerical  simulation  papers.  It  was  found  that  there  are 
three  methods  to  get  the  values  of  the  PEMFC  parameters.  One 
method  is  the  direct  measurement,  i.e.,  parameters  are  mea¬ 
sured  via  experiments.  For  example,  Lee  et  al.  [113]  prepared 
a  precise  impedance  measurement  system  based  on  two-probe 
and  four-probe  methods  to  measure  the  impedance  and  con¬ 
sequent  proton  conductivity  of  the  Nation  membrane.  Saito  et 
al.  [114]  measured  the  ionic  conductivity,  water  transference 
coefficient,  water  permeability  and  diffusion  coefficients  of  the 
water  and  the  Li+  cation  for  several  membranes.  Parthasarathy 
et  al.  [115]  and  Zhang  et  al.  [116]  measured  the  kinetic  and 
mass  transport  properties  for  the  oxygen  reduction  reaction  in 


the  membrane.  However,  some  parameters  are  very  difficult  to 
measure,  or  even  cannot  be  determined  by  experimental  meth¬ 
ods  and  must  be  estimated  [9].  Therefore  the  second  method 
to  obtain  the  parameters  is  via  an  appropriate  computational  or 
fitting  method.  Suares  and  Hoo  [117]  estimated  four  parame¬ 
ters  such  as  the  exchange  current  density  for  oxygen  reaction, 
diffusion  coefficient  of  water,  evaporation  and  condensation 
rate  and  overall  heat-transfer  coefficient  using  voltage-current 
data.  Carnes  and  Djilali  [80]  defined  an  algorithm  for  non¬ 
linear  least  squares  fitting  to  estimate  the  effective  membrane 
conductivity,  exchange  current  densities  and  oxygen  diffusion 
coefficients  in  a  one-dimensional  PEMFC  model.  Berg  et  al. 
[118]  estimated  four  parameters  such  as  exchange  current,  mem¬ 
brane  water  transfer  coefficient,  effective  oxygen  diffusivity  and 
average  membrane  resistance  using  a  one-dimensional  PEMFC 
model  based  on  [34].  Thamapan  et  al.  [119]  performed  a  param¬ 
eters  estimation  for  the  membrane  conductivity  submodels  using 
curve  fitting.  Guo  et  al.  [120]  fitted  cathode  catalyst  layer 
parameters  such  as  porosities,  reference  current  densities  and 
effective  diffusion  coefficients  using  a  one-dimensional  cath¬ 
ode  catalyst  layer  model.  However,  for  most  of  the  simulation 
researchers,  the  parameters  were  obtained  from  models  given 
by  other  research  groups,  which  constituted  the  third  method 
to  obtain  the  parameters.  Thus  it  can  be  seen  that  one  of  the 
major  reason  that  some  transport  and  electrochemical  param¬ 
eter  variation  ranges  are  quite  large  is  simply  because  of  the 
inherent  difficulty  and  complexity  in  experimental  measure¬ 
ment.  Actually  because  of  the  small  size  of  the  gas  channels  and 
the  small  thickness  of  the  diffusion  layers,  catalyst  layers  and 
membrane,  it  is  very  difficult  to  measure  the  flow  and  species 
distributions  in  these  regions.  This  situation  should  be  taken 
into  account  when  the  fuel  cell  model  validation  issue  is  con¬ 
cerned. 

From  previous  studies,  the  transport  and  electrochemical 
parameters  can  be  classified  according  to  their  influence  on  the 
model  results  as  insensitive  and  sensitive.  The  permeability  and 
solid  phase  conductivity  consist  of  the  insensitive  parameters. 
Most  of  the  other  parameters  such  as  diffusion  layer  porosity, 
membrane  phase  conductivity,  cathode  reference  exchange  den¬ 
sity  multiplied  by  area,  oxygen  diffusivity,  reference  oxygen 
concentration  et  al.  are  the  sensitive  parameters. 

As  far  as  the  validation  criterion  is  concerned  only  very 
recently  researchers  in  the  international  community  have  shown 
their  concern  on  whether  the  polarization  curve  only  is  enough 
for  the  model  validation  (the  CFD/NHT  model  validation).  Hak- 
enjos  et  al.  [121]  pointed  out  that  for  the  validation  of  mul¬ 
tidimensional  models,  using  only  the  polarization  curve  is  not 
sufficient,  and  they  performed  an  additional  comparison  between 
the  measured  and  simulated  electrical  current  distributions.  Lum 
and  McGuirk  [33]  also  used  a  two-step  validation  approach: 
global  validation  by  the  polarization  curve  and  local  validation 
by  the  distribution  of  local  current  density  obtained  from  a  seg¬ 
mented  fuel  cell.  In  [122]  an  interesting  example  was  presented: 
a  three-dimensional  PEMFC  model  was  used  for  a  single  chan¬ 
nel  fuel  cell.  In  one  case  the  ionic  resistance  in  two  catalyst 
layers  was  included,  while  in  the  other  case  these  resistances 
were  neglected.  By  adjusting  the  kinetics,  the  numerical  simu- 
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lation  results  for  the  current  density  of  the  two  cases  were  exactly 
the  same  at  a  voltage  of  0.75  V.  Thus  the  authors  proposed  that 
apart  from  the  global  validation,  the  local  distribution  of  current 
density  should  be  added  in  order  to  validate  a  comprehensive 
PEMFC  model. 

It  can  be  seen  that  after  about  20  years,  in  numerical  mod¬ 
eling  of  PEMFC  performance  researchers  have  become  aware 
that  the  polarization  curve  comparison  is  enough  for  the  vali¬ 
dation  of  a  comprehensive  model.  This  work  tries  first  to  per¬ 
form  a  comprehensive  study  of  the  parameter  sensitivity  for 
an  advanced  PEMFC  model  based  on  some  numerical  results 
to  give  a  detailed  discussion  of  how  to  validate  a  simulation 
model.  We  finally  found  that  even  the  proposed  two-step  vali¬ 
dation  approach  is  not  enough  to  validate  a  model,  hence  a  third 
validation  index,  which  is  relatively  easy  to  be  measured  is  pro¬ 
posed.  The  sensitivity  examination  results  and  the  discussion  of 
the  model  validation  issue  will  be  presented  in  the  companion 
paper  [123]. 

3.  Present  model  description 

The  computational  domain  of  the  present  model  is  shown 
in  Fig.  1.  The  conventional  parallel  flow  fields  are  adopted  in 
this  model.  The  model  assumes  that  the  fuel  cell  structure  is 
repeated  periodically  along  the  y-direction.  Neglecting  the  end 
effects  for  each  gas  channel,  it  can  be  regarded  that  the  pro¬ 
cess  in  each  channel  is  identical.  Hence,  to  save  computation 
time,  half  of  a  gas  channel  can  be  taken  as  the  computational 
domain  as  shown  in  Fig.  2.  Dry  air  is  fed  into  the  cathode 
channel,  whereas  humidified  hydrogen  is  supplied  to  the  anode 
channel. 

The  assumptions  adopted  in  the  present  model  are: 
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Fig.  2.  The  two-dimensional  cross-sections  of  the  computational  domain:  (a) 
y-z  cross-section;  (b)  x-z  cross-section. 


(1)  The  fuel  cell  operates  under  a  steady-state  condition. 

(2)  The  gas  mixture  is  an  incompressible  ideal  fluid. 

(3)  The  flow  in  the  gas  channels  is  laminar. 

(4)  The  diffusion  layer,  catalyst  layer  and  membrane  are 
isotropic  and  homogeneous,  and  the  membrane  is  consid¬ 
ered  impervious  for  reactant  gases. 

(5)  Ohmic  heating  in  the  bipolar  plates  and  the  diffusion  layers 
are  neglected  due  to  their  high  conductivities. 

(6)  Ohmic  potential  drops  in  the  diffusion  layers  and  bipolar 
plates  are  neglected  due  to  their  high  electrical  conductivi¬ 
ties. 

(7)  The  contact  resistance  between  any  two  parts  in  the  fuel  cell 
is  neglected. 


3.1.  Model  equations 


The  three-dimensional,  two-phase,  non-isothermal  model 
consists  of  non-linear,  coupled  partial  differential  equations  rep¬ 
resenting  the  conservation  of  mass,  momentum,  species,  charge 
and  energy.  The  conservation  equations  are  described  in  the  vec¬ 
tor  form  as  follows. 

Mass  conservation  equation: 


V  •  OgUg)  —  Sm 


Momentum  conservation  equation: 


1 


s(l  -  s ) 


V  •  (PgUgUg)  =  -V/?g  + 


1 


s(l  -  s ) 


OgVug)  +  Su 

(2) 


Species  conservation  equation: 

V  •  OgUgXk)  =  V  •  (pgDk,effVXk)  +  5k  (3) 

where  the  index  refers  to  different  species,  including  oxygen, 
hydrogen  and  water  vapor. 

Electrical  charge  equations: 


V  •  (/csV0s)  +  s  —  0 

(4) 

V  •  (/CmV0m)  +  —  0 

(5) 

Energy  conservation  equation: 

V  •  (pgUgT)  =  V  •  (AeffVT)  +  Sj 

(6) 

Source  terms  in  the  above  governing  equations  (5m, 

5u?  5k,  50?s, 

50,m,  and  Sj)  are  summarized  in  Table  2  for  various  sub-regions 
of  the  fuel  cell.  The  source  term  in  momentum  conservation 
equation,  5U,  represents  Darcy’s  drag  force  imposed  by  the  pore 
walls  on  the  fluid  within  the  pores,  which  usually  results  in  a 
significant  pressure  drop  across  the  porous  medium.  Eq.  (2)  is 
the  general  expression  of  the  momentum  equation.  In  the  gas 
channel  region,  the  porosity  e  becomes  unity  and  the  coefficient 
of  permeability  approaches  infinity,  hence  Eq.  (2)  resumes  a 
conventional  form  of  the  momentum  equation.  In  the  porous 
medium  region,  the  general  momentum  conservation  equation 
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Table  2 

Source  terms  for  governing  equations  in  various  regions  of  a  PEMFC 
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reduces  to  the  expression  of  Darcy’s  law: 

KKV. 

Ug  = - AVpg  (7) 

Jig 

KKri 

= - ~Vp\  (8) 

Ml 

The  source  term  in  mass  conservation  (5m)  and  species  conser¬ 
vation  equations  for  O2  and  H2  (5k)  are  the  volumetric  sink  or 
source  terms  due  to  the  electrochemical  reactions  in  the  cata¬ 
lyst  layer,  and  they  are  zero  in  other  parts  of  the  computational 
domain.  The  source  term  in  the  energy  conservation  equation 
(5t)  represents  the  sum  of  the  reversible  heat  release  and  irre¬ 
versible  heat  generation  [56]. 

In  Eqs.  (1)— (8)  and  in  Table  2  a  lot  of  parameters  and  variables 
need  to  be  further  determined.  They  are  described  as  follows. 

The  liquid  water  saturation  that  appears  in  Eq.  (2)  is  defined 
as  the  volume  fraction  of  liquid  water  in  the  porous  media,  that 
is 


V\ 

l -VS 


In  order  to  derive  the  governing  equation  for  the  liquid  water  sat¬ 
uration  the  mass  conservation  equation  of  liquid  water  is  needed: 


V  •  (yOUl)  =  -5W 


(10) 


The  so-called  capillary  pressure pc  is  defined  as  the  difference 
of  the  pressure  between  the  gas  and  the  liquid: 

Pc  =  Pg~P\  (Ha) 

This  pressure  is  assumed  to  be  a  function  of  liquid  water  satu¬ 
ration  [62,124]: 

pc  =  cf(55  [1.417(1  -  .?)  -  2.120(1  -  sf  +  1.263(1  -  sf] 

(lib) 


With  Eqs.  (1),  (7),  (8),  (10),  (11a)  and  (lib),  the  liquid  water 
saturation  equation  can  be  derived,  which  says: 


same  for  the  liquid  water  saturation  equation,  Eq.  (12),  and  it 
represents  the  interfacial  mass-transfer  rate  of  water  between  the 
gas  and  liquid  phases.  It  is  substantially  similar  to  the  form  in 
[39]: 


5W  —  hmss(pgXwsai  pgXw)  (13) 

where  hm  is  the  evaporation  and  condensation  rate  and  Xwsat 
is  the  mass  fraction  of  water  vapor  when  the  mixture  is  satu¬ 
rated  and  is  related  to  the  saturation  pressure  psat  at  operating 
temperature,  which  is  given  by  [36]: 

log10  psat  =  -2.1794  +  0.025937  -  9.1837  x  10"572 

+  1.4454  x  10“7r3  (14) 


The  determinations  of  a  number  of  diffusivities  are  now 
described.  The  diffusivity  in  Eq.  (3)  can  be  determined  as  fol¬ 
lows.  The  value  of  the  species  in  the  gas  channel  is  a  function 
of  temperature  and  pressure,  and  is  determined  by  following 
equation  [56]: 

7^k,eff  —  f^k,ref(^7  ^ref)  ^  (Pref/ P)  (13) 

where  Z\,ref  is  the  reference  value  at  Tref  and  pVQ f .  In  the  porous 
media  region  the  diffusivity  of  the  species  can  be  described  by 
the  Bruggeman  model  [56]: 

Ac,eff  =  £L5Dk,ref  (16) 


The  capillary  diffusion  coefficient  Dc  in  Eq.  (12)  is  given  by 
[82]: 


d pc 

Dc  = - ^ - — 

rj\  ds 


(17) 
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where  v  and  a  are  kinetic  viscosity  and  surface  tension,  respec¬ 
tively.  The  relative  permeabilities  for  the  liquid  and  gas  phases 
are  represented  by  [40] : 


V-  (ffl  — ^“gj  =  V  •  (Pi£>cW>  -  Sw  (12) 

V  m  KIg  V 

The  source  term  for  the  water  vapor  equation  (included  in 
Eq.  (3)  when  the  index  k  takes  the  corresponding  value)  is  the 


Kd  =  s 3,  £rg  =  (l-s)3  (19) 

Now  attention  is  turned  to  the  electrical  charge  equations.  The 
source  terms  in  Eqs.  (4)  and  (5)  are  directly  related  to  the  elec¬ 
trochemical  reaction  expressed  by  the  electrical  current,  which 
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is  given  by  Bulter-Volumer  equation: 
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(21) 


The  proton  conductivity  km  in  Eq.  (5)  is  related  with  the  water 
content  of  the  membrane,  X,  which  is  in  turn  a  function  of  the 
water  activity,  a  [20] : 


A"m 


(0.5139A. 


0.326)exp 


(22) 


0.043  +  17.81a  —  39.85a2  +  36.0a3  0  <  a  <  1 

14  +  1.4(a  —  1)  1  <  a  <  3 

(23) 


a  = 


cowRT 


pSSLt 


(24) 


where  cow  is  the  molar  fraction  of  water  vapor. 

Once  the  membrane  phase  potential,  0m  and  the  proton  con¬ 
ductivity  on  the  membrane,  km  are  obtained,  local  current  den¬ 
sity,  /,  can  be  calculated  by 


—  Adn  V  0m 

(25) 

The  overpotential  is  described  as 

=  Otot  10s  0s, ref  1  10m  0m, ref  1 

(26) 

where  rjtot  is  the  total  overpotential  of  anode  or  cathode,  0s?ref 
the  solid  phase  potential  at  reference  state  and  0m,ref  is  the  mem¬ 
brane  phase  potential  at  reference  state.  For  the  solid  phase 
potential,  the  potential  at  the  interface  between  the  anode  current 
collector  and  the  diffusion  layer  is  zero,  and  for  the  membrane 
phase  potential,  the  potential  at  the  interface  between  the  anode 
catalyst  layer  and  the  membrane  is  set  to  be  zero  also. 

The  operating  potential  of  the  cell  is  then  calculated  by 


Tcell  —  Voc  7?a,tot  7?c,tot  ^?m,pro 


(27) 


where  ^m,pro  is  the  Ohmic  overpotential  in  the  membrane,  and 
yoc  is  the  open  circuit  potential,  which  is  calculated  by  [21]: 


Voc  =  1.23-0.9  X  10_3(r 


298) +  2.3— log  (pipe) 

(28) 


Table  3 


Grid-independence  test 


Grid  size 

4v  (A  cm  2) 

12  x  12  x  40 

0.8034 

22  x  12  x  40 

0.8045 

32  x  12  x  40 

0.8058 

32  x  22  x  40 

0.8072 

42  x  22  x  40 

0.8074 

32  x  32  x  40 

0.8080 

It  should  be  noted  that  although  the  model  equations  of 
the  present  model  are  mainly  copied  from  existing  papers,  the 
present  two-phase  model  is  different  from  all  of  the  existing 
models.  The  previous  papers  known  to  the  authors  presented 
either  an  isothermal  one-dimensional  model  [20],  or  isothermal 
two-dimensional  models  [36,39,40,62,124],  or  a  single-phase 
three-dimensional  model  [56]  or  an  isothermal  two-phase  three- 


Fig.  3.  Flow  diagram  of  the  solution  procedure. 
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Fig.  4.  Polarization  curve  of  PEMFC  on  the  base  case:  (a)  Veen  =0.8;  (b) 
Ken  =0.3. 
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dimensional  model  [82].  However,  the  present  model  is  a  three- 
dimensional,  two-phase,  non-isothermal  model.  Furthermore, 
the  present  model  presents  a  method  in  detail  to  obtain  the 
voltage  versus  current  curve  of  PEMFC.  This  method  will  be 
described  in  Section  3.3. 

3.2.  Boundary  conditions 

In  the  x-z  plane,  symmetrical  conditions  are  adopted.  That 
is,  the  gradient  in  the  y-direction  of  each  variable  is  zero. 

At  the  gas  channel  inlet,  the  temperature  and  gas  species 
concentrations  are  assumed  to  be  uniform.  The  inlet  velocities 
are  specified  by 


u 


a,m 


I rpf 

=  U  —  A 


RT 


a,m 


1  1 


U 


c,m 


=  ?< 


IF 

Aef 


m 


Pa,i 


m 


A 


RT( 


c,m 


o,m 

1 


Ach 

1 


4F  Pc, in  2fh,in  Ach 


(29) 


(30) 


where  fa  and  are  the  reactant  stoichiometric  flow  ratio  of 
anode  and  cathode,  respectively,  they  are  defined  as  the  ratio 
of  the  amount  of  reactant  supplied  to  the  amount  of  reaction  to 
generate  the  specified  reference  current  density  /ref.  Am  is  the 
geometrical  area  of  the  membrane  and  Ac h  is  the  cross-sectional 
area  of  the  gas  channel. 
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Fig.  5.  Oxygen  mass  fraction  distribution  in  the  cathode:  (a)  Vceii=0.8;  (b) 
Kell  =0.3. 


Fig.  6.  Local  current  density  distribution  in  the  cathode  catalyst  layer:  (a) 
Ken  =  0.8;  (b)  Ken  =0.3. 
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A  local  one-way  assumption  is  adopted  to  give  the  gas  chan¬ 
nel  outlet  velocity  condition,  and  is  then  corrected  by  a  global 
mass  conservation  constraint  [125]. 

At  the  body  surface,  the  no-slip  condition  is  applied  for  the 
velocity  and  non-permeable  condition  of  the  species  mass  frac¬ 
tion. 

For  the  liquid  water  saturation,  the  computational  domain 
involves  two  diffusion  layers,  two  catalyst  layers  and  a  mem¬ 
brane.  At  the  interface  between  the  gas  channel  and  the  diffusion 
layer,  the  liquid  water  velocity  is  set  to  zero. 

The  computational  domain  for  the  electrical  charge  equations 
involves  the  anode  catalyst  layer,  the  membrane  and  the  cathode 
catalyst  layer.  The  boundary  conditions  are  described  as  follows: 

at  the  surface  of  the  anode  catalyst  layer  :  0S  =  0, 

90m  _  q 

dz 

at  the  surface  of  the  cathode  catalyst  layer  :  0S  =  Vceii, 


dz 


3.3.  Numerical  procedures 

The  governing  equations,  together  with  the  boundary  condi¬ 
tions  are  discretized  by  the  finite  volume  method.  The  SIMPLEC 
algorithm  [125]  is  utilized  to  deal  with  the  coupling  of  the  veloc¬ 
ity  and  the  pressure.  Since  all  governing  equations  are  coupled 
with  each  other,  they  ought  to  be  solved  simultaneously  with 
an  iterative  method.  The  solution  is  considered  to  be  convergent 
when  the  relative  error  of  each  dependent  variable  between  two 
consecutive  iterations  is  less  than  1.0  x  10-5. 

The  grid  system  used  is  32  x  22  x  40.  To  simulate  the 
local  transport  phenomena  in  the  fuel  cell,  the  grid  arrange¬ 
ment  at  z-direction  is  non-uniform.  The  grid-independence  test 
is  performed  on  six  grid  systems.  The  results  of  the  aver¬ 
age  current  density  computed  by  the  model  under  different 
grid  systems  when  the  fuel  cell  operating  voltage  is  0.5  V  are 
summarized  in  Table  3.  Considering  both  accuracy  and  eco¬ 
nomics,  the  grid  system  of  32  x  22  x  40  was  selected  for  present 
research. 

There  are  generally  two  ways  of  obtaining  the  voltage  versus 
current  curve:  either  the  operating  current  density  is  given 
and  different  potential  losses  are  calculated,  or  the  so-called 


Table  4 


Model  parameters  for  basic  case 


Parameter 

Symbol 

Value 

Reference 

Gas  channel  length 

L 

0.04  m 

Gas  channel  width 

W 

7.62  x  10-4  m 

Gas  channel  height 

Hch 

7.62  x  10-4  m 

Diffusion  layer  height 

Hd 

2.54  x  10“4  m 

Catalyst  layer  height 

2.87  x  10~5  m 

Membrane  height 

Hm 

2.3  x  10~4  m 

Land  area  width 

wcc 

7.62  x  10-4  m 

Faraday’s  constant 

F 

96487  C  mor1 

Gas  channel  inlet  temperature 

Tin 

353  K 

Anode/cathode  pressure 

PdPc 

1/1  atm 

Electron  number  of  anode  reaction 

na 

4 

Electron  number  of  cathode  reaction 

nc 

2 

Fuel/air  stoichiometric  flow  ratio 

£a/£c 

3/3 

[48] 

Relative  humidity  of  inlet  fuel 

RHa 

100% 

[56] 

Relative  humidity  of  inlet  air 

RHC 

0 

[56] 

Oxygen  mass  fraction  of  inlet  air 

0.23 

H2  diffusion  coefficient  at  reference  state 

Dh.ref 

0.915  x  10-4  m2  s_1 

[48] 

O2  diffusion  coefficient  at  reference  state 

f^o,ref 

0.22  x  10-4  m2  s_1 

[48] 

Water  vapor  diffusion  coefficient  at  reference  state 

f^w,ref 

0.256  x  10-4  m2  s_1 

[48] 

Anode  exchange  current  density  multiply  specific  area 

Asia, ref 

5.0  x  107  Am-3 

Cathode  exchange  current  density  multiply  specific  area 

A  s^'c, ref 

120  Am-3 

Hydrogen  reference  concentration 

Ch,ref 

56.4  mol  m-3 

[21] 

Oxygen  reference  concentration 

Co, ref 

3.39  mol  m-3 

[21] 

Anode  transfer  coefficient 

Oia 

0.5 

[56] 

Cathode  transfer  coefficient 

ac 

0.5 

[56] 

Porosity  of  diffusion  layer 

£d 

0.3 

[68] 

Porosity  of  catalyst  layer 

£ct 

0.28 

[81] 

Absolute  permeability 

K 

1.76  x  10“n  m2 

[36] 

Solid  phase  conductivity 

Ks 

53  Sm'1 

[V] 

Membrane  phase  conductivity 

Km 

6Sm_1 

[48] 

Surface  tension 

a 

0.0625  N  m-1 

[40] 

Evaporation  and  condensation  rate 

h  m 

7 

t/3 

O 

O 

H 

[80] 

Current  collector  thermal  conductivity 

kc 

150  Wm_1  K-1 

[36] 

Diffusion  layer  thermal  conductivity 

kd 

150  Wm_1  K-1 

[36] 

Membrane  thermal  conductivity 

km 

0.95  Wm_1  K-1 

[56] 

W.Q.  Tao  et  al.  /  Journal  of  Power  Sources  160  (2006)  359-373 


369 


potentiostatic  approach  is  used,  where  the  cell  potential  is  set 
and  the  current  density  is  calculated  [36].  We  chose  the  second 
approach  for  the  simulation.  By  giving  the  initial  values  of 
the  anode  total  overpotential  and  cathode  total  overpotential, 
the  current  density  can  be  obtained,  and  the  Ohmic  losses 
can  be  calculated.  Then  the  anode  total  overpotential  and 
cathode  total  overpotential  are  corrected  by  the  constraint  that 
the  anode  current  should  be  equal  to  cathode  current,  that 
is 


y^(4Vcv) 


YflcVcv) 


(31) 


where  /av  and  Vcv  are  the  cell  average  current  density  and  the 
volume  of  the  control  volume,  respectively.  Such  an  iterative 
solution  procedure  is  shown  in  Fig.  3. 

To  conduct  the  modeling  simulation,  a  great  number  of 
parameters  are  required.  The  parameters  of  the  basic  case  and  the 
corresponding  references  are  all  listed  in  Table  4.  The  numerical 
results  of  the  basic  case  are  regarded  as  the  references  for  further 
parameter  sensitivity  examination. 


4.  Results  and  discussion  for  the  basic  case 

In  this  section,  the  PEMFC  polarization  curve  for  the  basic 
case  will  be  first  presented.  Then  the  distribution  of  the  oxygen 
mass  fraction  in  the  cathode,  local  current  density,  liquid  water 
saturation  and  temperature  will  be  presented  in  order. 

The  polarization  curve  of  the  basic  case  is  shown  in  Fig.  4.  It 
follows  the  general  variation  trend  observed  for  the  PEMFC  in 
many  references  and  will  not  be  further  discussed. 

Fig.  5  shows  profiles  for  the  oxygen  mass  fraction  in  the 
cathode,  including  the  gas  channel  and  diffusion  layer.  At  high 
cell  voltages,  the  oxygen  mass  fraction  is  relatively  uniform. 
Whereas  at  low  cell  voltages,  the  oxygen  mass  fraction  is  far 
from  being  uniform,  which  implies  that  the  local  current  den¬ 
sity  is  non-uniformly  distributed  in  the  catalyst  layer  at  a  low 
cell  voltage  since  the  local  current  density  is  dependent  on  the 
oxygen  concentration. 

Fig.  6  shows  the  distribution  of  local  current  density  in  the 
cathodic  catalyst  layer.  It  can  be  seen  that  the  distribution  is 
quite  uniform  at  high  cell  voltages.  The  local  current  density  is 
somewhat  lower  over  the  shoulder  than  over  the  channel.  The 
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Fig.  7.  Liquid  water  saturation  distribution  in  the  cathode  electrode  (Vceii  =  0.6): 
(a)  x-y  plane;  (b)  y-z  plane. 


Fig.  8.  Liquid  water  saturation  distribution  in  the  cathode  electrode  (VCeii  =  0.3): 
(a)  x-y  plane;  (b)  y-z  plane. 
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minimum  current  density  is  located  at  the  corner  of  the  catalyst 
layer  over  the  shoulder  adjacent  to  the  diffusion  layer.  While  at 
low  cell  voltages,  the  distribution  pattern  becomes  very  different 
from  that  at  high  cell  voltages.  The  minimum  current  density 
is  located  at  the  corner  of  the  catalyst  layer  over  the  shoulder 
adjacent  to  membrane,  and  the  current  density  distribution  is 
very  non-uniform.  The  phenomena  result  from  the  mass  transfer 
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Fig.  10.  Cathodic  and  anodic  overpotentilas  vs.  cell  voltage. 

limitation  of  oxygen.  The  oxygen  concentration  decreases  along 
the  z-direction  close  to  the  membrane  due  to  consumption  and 
mass  transport  resistance. 

Now  the  liquid  water  saturation  in  the  cathode  diffusion 
layer  and  catalyst  layer  is  discussed.  The  results  are  shown  in 
Figs.  7  and  8.  Fig.  7  shows  that  at  high  cell  voltages,  the  liquid 
water  saturation  is  lower  and  relatively  uniform,  and  it  increases 
gradually  along  the  direction  from  the  gas  channel  inlet  to  the 
outlet.  The  maximum  liquid  saturation  is  located  at  the  corner  of 
the  catalyst  layer  over  the  shoulder  adjacent  to  the  membrane. 
Fig.  8  shows  that  at  low  cell  voltages,  the  liquid  water  saturation 
distribution  is  very  different  from  that  at  high  cell  voltages.  The 
liquid  water  saturation  in  the  domain  vary  significantly,  and  the 
maximum  value  is  not  at  the  corner  of  the  catalyst  layer  over  the 
shoulder  adjacent  to  membrane  but  at  the  corner  of  catalyst  near 
the  gas  channel  inlet.  The  reason  is  that  at  that  point,  the  oxygen 
mass  fraction  reaches  a  maximum  value  and  the  generated  liquid 
water  cannot  be  extracted  in  time. 

Fig.  9  shows  the  temperature  distribution  in  the  PEMFC  when 
the  cell  voltage  is  0.6  V.  The  results  show  that  the  temperature  at 
the  cathode  side  is  higher  that  at  the  anode  side.  The  maximum 
temperature  located  at  the  cathode  catalyst  layer  over  the  gas 
channel  is  due  to  reversible  and  irreversible  entropy  production. 
The  maximum  temperature  increase  is  about  6  °C.  Beming  et 
al.  [48]  simulated  the  maximum  temperature  increase  as  about 
3  °C  when  the  average  current  density  is  1.2  A  cm-2.  Whereas 
the  result  modeled  by  Ju  et  al.  [56]  was  over  10  °C  at  Vceii  =  0.6. 
This  paper  gives  an  intermediate  result. 

The  variations  of  the  cathodic  and  anodic  overpotentials  with 
the  current  density  are  presented  in  Fig.  10.  It  is  interesting  to 
note  that  in  most  of  the  previous  studies,  such  information  was 
usually  not  provided.  To  the  authors’  knowledge,  only  Ref.  [38] 
presented  such  information.  Our  results  show  [123]  that  such 
information  is  useful  for  the  validation  of  a  model. 

5.  Conclusions 


Fig.  9.  Temperature  distribution  in  the  PEMFC  (Vceii  =  0.6):  (a)  x—y  plane;  (b)  A  comprehensive  review  of  PEMFC  models  in  the  Open  lit- 

y-z  plane;  (c)  x-z  plane.  erature  was  conducted.  From  this  review  it  was  found  that  at 
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least  10  different  models  cited  the  same  test  data  for  the  V-I 
curve  to  verify  their  correctness,  while  the  physical  and/or  elec¬ 
trochemical  parameters  involved  in  different  models  were  quite 
different,  not  only  different  in  quantity  but  also  different  in  the 
order  of  magnitude.  This  situation  stimulated  the  present  authors 
to  perform  a  sensitivity  study  for  the  major  parameters  and  to 
examine  whether  the  V-I  curve  only  can  be  serve  as  the  model 
verification  index. 

A  three-dimensional,  two-phase  and  non-isothermal  model 
was  developed  based  on  the  existing  models.  The  simulated 
results  include  the  polarization  curve,  the  oxygen  mass  fraction 
distribution  in  the  cathode,  the  local  current  density  distribu¬ 
tion  in  the  cathode  catalyst  layer,  the  liquid  water  saturation 
distribution  in  the  cathode  electrode,  the  cathodic  and  anodic 
overpotentials  versus  current  density,  and  the  temperature  dis¬ 
tribution  in  the  PEMFC.  Generally  speaking,  these  simulated 
results  qualitatively  agree  with  the  existing  results  in  the  litera¬ 
ture. 

The  parameter  sensitivity  examination  results  and  discussion 
of  model  validation  are  reported  in  the  companion  paper. 
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